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Abstract. A new two-dimensional axi-symmetric ultraviolet radiative transfer code is presented, which is used 
to calculate photodissociation and ionization rates for use in chemistry models of flaring circumstellar disks. 
Scattering and absorption of photons from the central star and from the interstellar radiation field are taken 
into account. The molecules are effectively photodissociated in the surface layer of the disk, but can exist in the 
intermediate, moderately warm layers. A comparison has been made with an approximate 2D ray-tracing method 
and it was found that the latter underestimates the ultraviolet field and thus the molecular photodissociation 
rates below the disk surface. The full 2D results show significantly higher abundances of radicals such as CN and 
C 2 H than previous work, partly due to the fact that CO is dissociated to greater depths. Results for different 
stellar radiation fields are also presented. The CN/HCN ratio shows a strong dependence on the stellar spectrum, 
whereas other ratios such as HCO + /CO show only little variation. 
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1. Introduction 


Numerous recent observations of dust and molecular line 
emission have increased the understanding of disks around 
isolated T-Tauri stars in both their appearances and in¬ 
ternal structure (Beckwith & Sargent 1996, Dutrey et 
al. 2000, Qi 2001). The modeling to date can be di¬ 
vided in two categories. A first category deals with 
the physical structure of disks. It was proposed by 


Kenyon fe Hartmann (1987) that disks have a flared ge¬ 
ometry, where the stellar light puffs up the outer regions 
of the disk. Recent papers by Chiang fe Goldreich (1997, 
1999), Men’shchikov & Henning (1997), D’Alessio et al. 


(1998, 1999), Dullemond et al. (2002) and Nomura (2002) 
show similar results by calculating explicitly the density 
and temperature distributions in disks in hydrostatic equi¬ 
librium. Chiang & Goldreich (1997) fit spectral energy 
distributions by considering disks with a bi-layered tem¬ 
perature structure. This structure is calculated assuming 
that the ultraviolet and visible light from the star heats 
up the outer layer of the disk, and does not effect the 
shielded region directly. D’Alessio et al. (1998) calculate 
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the full temperature structure in the vertical direction and 
include the effects of accretion and viscous heating. 

The second category takes a physical disk model and 
uses a time-dependent chemistry code to calculate the 
molecular abundances in the gas (Aikawa et al. 1997; 
Aikawa & Herbst 1999; Willacy et al. 1998; Willacy & 
Langer 2000). Aikawa & Herbst (1999) use the minimum- 
mass solar nebula model of Hayashi (1981), whereas 
Willacy et al. (2000) use the Chiang & Goldreich model. In 
the first model, the temperatures are very low throughout 
the disk, thereby freezing out the molecules onto grains, 
while in the latter model nearly all molecules are disso¬ 
ciated due to the ultraviolet light in the upper layer and 
are depleted in the cold lower layer. Both models have to 
assume either an artificially low sticking probability or a 
very high photo-desorption rate to keep molecules in the 
gas. 

In the latest paper by Aikawa et al. (2002, hereafter 
Paper I), the D’Alessio et al. model was used to calculate 
molecular abundances. The stellar and interstellar ultra¬ 
violet (UV) radiation heats the dust and gas, and dis¬ 
sociates molecules in the surface region of A v <; 1 mag. 
In the deeper layers, visible photons with lower energies 
are available to heat the grains and gas, but there are 
much fewer ultraviolet photons for dissociation than in 






2 


van Zadelhoff et al.: Chemical effects due to 2D UV transfer 


the surface layer. The dust in the surface layer, heated by 
the ultraviolet and visible photons, in turn emits infrared 
photons, which can travel even further toward the mid¬ 
plane until they are eventually absorbed. Viscous heating 
is only important in the inner few AU of the disk, which 
are not considered in this study. As a result, the disk has 
three distinctive chemical layers. In the surface layer with 
temperatures T > 40 K and density un ~ 10 4 — 10 5 
cm -3 , molecules are mostly dissociated. In the interme¬ 
diate layer, conditions are comparable to those in dense 
interstellar clouds, i.e., the temperature is 20-40 K and 
density is ~ 10 6 — 10' cm" 3 . This layer contains sufficient 
amounts of gaseous molecules to reproduce the observa¬ 
tions of some disks. In the innermost layer including the 
midplane, the density is high (% > 10 7 cm -3 ) and tem¬ 
perature is low (T <; 20 K). Nearly all the gaseous species, 
except for H 2 and He, freeze out onto the grains and stay 
there until the grains are again heated, e.g. when most 
of the disk has dissipated and becomes optically thin to 
ultraviolet radiation or when the grains are dynamically 
transported to a warmer part of the disk. 

In this work, we again adopt the density and temper¬ 
ature distribution of the D’Alessio et al. model, and use 
a new two-dimensional (2D) ultraviolet radiative trans¬ 
fer code to describe the radiative effects on the chemistry. 
Previous works, including our Paper I, treated the UV and 
optical radiative transfer in disks in a one dimensional 
(ID) or approximate 2D method. In the former case it 
is assumed that the disks are geometrically thin, i.e. the 
scale-height of the disk is small compared to the distance 
to the star, and therefore radiative transfer in the radial 
direction is negligible compared to that in the vertical di¬ 
rection (Chiang & Goldreich 1997, D’Alessio et al. 1998). 
In the latter case the attenuation both along the line of 
sight from the star and in the vertical direction is cal¬ 
culated (Paper I). This method cannot treat scattering 
effects in principle and therefore only gives an approx¬ 
imate depth-dependent intensity; forward scattering will 
help stellar photons at UV and optical wavelengths to pen¬ 
etrate deeper into the disk, which will affect the chemistry 
through photodissociation. One important aim of this pa¬ 
per is a direct assessment of the effects of the correct 2D 
treatment compared with the earlier more approximate 
methods. The 2D radiative transfer code also enables us 
to consider dissociation of H 2 and CO via the stellar UV 
radiation, which was neglected in Paper I. In addition the 
dependence of molecular abundances and line intensities 
on the stellar spectrum can be studied. 

Paper I and this paper are complementary to the re¬ 
cent work of Markwick et al. (2002). Whereas our studies 
focus on the outer regions of the disk (>50 AU), Markwick 
et al. (2002) consider the inner 10 AU in a quasi-static 
treatment; the chemistry is calculated locally, with the 
starting conditions defined by the chemical abundances 
calculated at more distant radial points. Their calcula¬ 
tions confirm the results of earlier ID studies that the 
chemistry in the inner region is driven by thermal desorp¬ 
tion. The inclusion of X-rays in their model shows that 


radicals such as CN and C 2 H are important tracers of the 
ionization in disks, a conclusion also reached for the outer 
disks in our papers. 

The outline of this paper is as follows. In the cou¬ 
pling between the radiation and chemistry is discussed. In 
the results are presented of the inclusion of continuum 
radiative transfer on the chemistry and compared with the 
approximate 2D method. The differences in line strength 
of the various molecules are calculated in §||. Discussion is 
in §||, followed by a conclusion in §6. Finally the adopted 
2D Monte Carlo UV radiative transfer code is explained 
in detail in Appendix |a|. 


2. Chemistry and ultraviolet radiative transfer 

2.1. General considerations 

The gas-phase molecular abundances are governed by two 
different effects of the radiation. On the one hand, the stel¬ 
lar radiation heats up the disk, assuring that freeze-out 
onto grains does not occur. This heating by direct radia¬ 
tion is dominated by the lower energy photons emitted by 
the cool central star (T* = 4000 K). In this paper, we as¬ 
sume that the gas and dust temperatures are coupled and 
adopt the calculated dust temperatures from D’Alessio et 
al. (1999). On the other hand, the higher-energy UV ra¬ 
diation dissociates and/or ionizes molecular and atomic 
species and most T Tauri stars have a relatively large UV 
flux compared to main sequence stars of the same tem¬ 
perature (e.g. Costa et al. 2000). The mean intensity of 
the visible and ultraviolet light at a position in the disk 
is controlled by absorption and scattering, making the so¬ 
lution non-trivial. In this section, we briefly discuss the 
UV radiative transfer in the disk. The numerical details 
of the Monte Carlo code are explained in more detail in 
Appendix [\j. 

Before combining the radiative transfer with the chem¬ 
istry, the different time scales in the problem are checked. 
The three time scales of interest are: the radiative time 
scale, the chemical time scale, and the dynamical time 
scale. Since the current line observations of disks have poor 
spatial resolution, even for sources mapped with millime¬ 
ter interferometers (Dutrey et al. 1996, Qi 2001), they are 
mostly sensitive to the outer regions of the disk. A com¬ 
parison of the time scales is therefore most appropriate at 
distances > 100 AU from the star. 

The first time scale represents the time needed for the 
radiation field to adapt to fluctuations in physical and/or 
chemical processes. The disk itself is too cold to emit ul¬ 
traviolet or visible photons, therefore the mean intensity 
at each point in the disk depends only on changes in the 
source, both stellar or interstellar, and the density distri¬ 
bution in between the source and that point. The radiation 
field reaches its mean value on the same time scale as it 
takes the light to travel from the source to the point of in¬ 
terest; for 100 AU this is ~ 14 hours from the star. So the 
radiation field is always in equilibrium in the time scales 
we are concerned with in this paper. 
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The chemical time scale, i.e., the time in which the 
number density of a species reaches its steady-state value, 
varies with molecular species and position in the disk. In 
the intermediate molecular layer we are concerned with, it 
is of the order of 10 4 -10 5 years for many organic molecules 
(Paper I). On the other hand, in the midplane layer the 
time scale of adsorption, which is the dominant process, 
is very short: ~ 10 2 (10 8 cm -3 /n.H) yr. 

There are several dynamical time scales in the disk. 
The typical accretion time scale of the disk as a whole 
is estimated to be 10 6 yr, derived from dividing the disk 
mass (10“ 2 Mq) by the mass accretion rate (10 -8 M 0 
yr -1 ). The local accretion timescale at R = 100 AU, 
i.e. the time it takes for a parcel of gas to move 5 AU 
(less than the radial gridsize in our model) radially in¬ 
wards, is estimated to be shorter, 7 x 10 4 yr from the 
local radial drift velocity wr = M/[2ttRT,(R)\. Since ac¬ 
cretion is caused by the transport of angular momentum 
due to turbulence, the timescale for turbulent mixing is 
also relevant. This timescale for radial mixing over A R 
is ~ (Ai?) 2 /(/ t urbfturb), for A R larger than the size of 
turbulent eddy / t .urb (e.g. Aikawa et al. 1996). The actual 
Eddy size (Zturb) and turbulent velocity (i>t ur b) are highly 
uncertain, but if they are assumed to be 10% of the verti¬ 
cal scale height and sound velocity, the timescale for global 
mixing (i.e. A R ~ R ~ 100 AU) is of the order of 10 6 yr. 
The local mixing timescale with smaller A R (or A Z) is 
accordingly shorter. 

In this paper we assume that the disk is static; the 
molecular abundances are calculated at each point ( R , Z) 
in the disk, without considering any accretion flow or mix¬ 
ing for simplicity. Although this static model is invalid 
from a physical point of view, as the comparison of the 
different local time-scales show above, this simplification 
is adopted for the following two reasons. First, disk chem¬ 
istry is a complicated mixture of various chemical pro¬ 
cesses and physical processes (e.g., density and temper¬ 
ature distributions, radiative transfer and hydrodynam¬ 
ics). In this series of papers, the importance of each of 
the processes is investigated step by step; in Paper I the 
chemical processes which are important in the different re¬ 
gions of the disk were analyzed, and in the current paper 
the effect of 2-dinrensional UV radiative transfer is inves¬ 
tigated. Mixing and accretion smooth out the molecular 
stratification over some length, determined by the balance 
between the chemical and dynamical time scales. At the 
same time radical species are transfered in the interme¬ 
diate layer, activating chemical reactions in neighbouring 
regions. These effects are expected to be larger in the ver¬ 
tical direction than in the radial direction, since the radial 
distributions do not vary much in the outer disk between 
50 and 400 AU (Paper I). They should be investigated in 
detail in independent future work. Second, full coupling of 
these effects in 2- or 3-dimensional calculations would be 
very complicated and too time-consuming with current 
computational facilities. In order to accomplish the full 
coupling of chemistry, radiative transfer and (magneto- 
hydrodynamics, the chemical reaction network needs to 


be significantly reduced by carefully choosing the domi¬ 
nant chemical reactions in each region of the disk. The 
results presented in our papers using a large chemical net¬ 
work will assist the construction and testing of future cou¬ 
pled chemical-hydrodynamical models. 

Another relevant physical process is the settling of 
dust and grain-growth. In this paper we assume for sim¬ 
plicity that the optical properties of dust grains are the 
same as those of interstellar clouds, and that the gas and 
dust are well mixed. Settling of small (sub-micron sized) 
grains, which are the main source of the UV opacity, takes 
place over a period of at least 10 6 years (e.g. Nakagawa, 
Nakazawa, & Hayashi 1981; Weidenschilling 1997). The 
coagulation timescale, however, is much shorter. In fact, 
the spectral energy distributions (SEDs) of T Tauri stars 
indicate significant increase to mm-sized grains compared 
with interstellar dust (Chiang et al. 2001; D’Alessio, 
Calvet, & Hartmann 2001). Inclusion of dust coagulation 
in our radiation-chemistry disk model is however more 
complicated, since variations in dust-properties, scatter¬ 
ing and possible decoupling of the gas and dust tempera¬ 
tures could be of essential importance. In previous work on 
the vertical disk structure, the temperature distribution is 
determined by radiative transfer, with dust as the main 
source of opacity, whereas the density distribution is deter¬ 
mined by the hydrostatic equilibrium, assuming the same 
gas and dust temperatures. If the number of small grains 
decreases due to coagulation, the two temperatures may 
decouple if not enough collisions between gas and dust 
particles take place. Calculations of the decoupled gas and 
dust temperature have only been performed for tenuous 
disks which have dust masses that are at least two orders 
of magnitude less than the disks studied in this work and 
are optically thin to UV and optical radiation, simplifying 
the radiative transfer (Kamp & van Zadelhoff 2001). This 
problem of dust coagulation is left for future work, and 
a basic 2-D radiation-chemistry model with interstellar¬ 
sized grains is established in this paper. 

The chemical model adopted in this paper is the 
“new standard model” (NSM) for the gas-phase chemistry 
(Tervieza & Herbst 1998; Osamura et al. 1999), extended 
to include deuterium chemistry, adsorption of molecules 
onto grains, and thermal desorption (Aikawa & Herbst 
1999; and references therein). Chemical processes by X- 
rays are not included in this paper. 

2.2. Impact of ultraviolet radiation 

The chemical abundances are related through a reaction 
network. In this network the dissociation- and ionization- 
rates ki of the different species are calculated using 

f°° X 

ki — / 47 t — cri{\)J\d\ ( 1 ) 

Jo " c 

where h is the Planck constant, c the velocity of light and 
J\ the mean intensity at wavelength A, defined as 

Jx = ~~ f Ix dn. ( 2 ) 

4tt J n 
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Some molecules can be dissociated by absorption of pho¬ 
tons over a continuous part of the spectrum, while others 
need discrete line transitions or a combination of the two 
(van Dishoeck 1988). The cross sections <7* (A) or oscilla¬ 
tor strengths /((A) need to be specified for each transi¬ 
tion for each molecule. We adopt here the values given by 
van Dishoeck (1988), updated by Jansen et al. (1995a & 
1995b). For species for which no data on the cross sections 
are available, a rate given by a similar type of molecule 
was adopted. 

The photodissociation rate of a species depends on 
the mean intensity of the radiation in the spatial grid, 
which is obtained via the 2D Monte Carlo code. This 
on the same principles as those developed by 
, Spaans (1996) and Gordon et al. (2001) for 
interstellar applications and e.g., Whitney & Hartmann 
(1992) for circumstellar disks. In general, one needs to de¬ 
fine a source with a surface area S, where the total energy 
is the flux from the source times the surface area. The total 
energy is divided into a number of photon-packages (N\), 
which then travel through the grid. The energy density at 
each point is calculated using both the number of photons 
passing through a grid cell and the mean energy trans¬ 
ported by each photon during its transfer. Photon pack¬ 
ages are either absorbed or scattered by dust grains and 
are followed until they leave the system or have reached a 
minimum intensity. Optical properties of dust grains are 
assumed to be the same as those of interstellar dust; the 
scattering is described by the Henyey-Greenstein probabil¬ 
ity function (Eq. |A.8| ), and albedo and opacities are taken 
from Roberge et al. (1991). 

In the Monte Carlo method, the radiation field within 
the disk is readily obtained for irradiation by both stel¬ 
lar and interstellar radiation: since each source of radia¬ 
tion gives its own energy density, calculations for different 
sources can be added to give a total energy density. In 
addition, radiation fields for different stellar spectra are 
easily calculated (Sec. 3.4). The energy densities in the 
disk are directly proportional to the emitted stellar flux, 
and are thus scalable. Since the disk structure is indepen¬ 
dent on the precise UV spectrum of the central star, a 
single calculation is sufficient for obtaining the full UV 
radiation distribution for each assumed spectrum. 

In the 912-1110 A region, additional attenuation due 
to H 2 and CO self-shielding was included. The shield¬ 
ing of other species by H 2 was ignored in these calcula¬ 
tions. This could lead to an overestimate of dissociation 
rates for species which dissociate only in the 912-1110 A 
region by at least 20 % (Draine & Bertoldi 1996). Self¬ 
shielding of H 2 and CO depends on the column densities 
of the molecules themselves along the photon trajectories 
as well as the dust attenuation (A v ), which means that 
the chemical abundances and the UV radiative transfer 
need to be solved simultaneously. For the interstellar UV 
field this is easily done via a ID slab model (van Dishoeck 
& Black 1988, Lee et al. 1996), which is adopted in this 
paper and previous works (e.g. Paper I). Shielding of the 
stellar UV radiation is more difficult to calculate, since it¬ 


code is based 
Boisse (1990) 


eration of the chemistry and 2D radiative transfer is very 
time-consuming. Since we are not interested in the precise 
position of the H — H 2 transition, this iteration is per¬ 
formed only once: we first calculate the mean intensity as¬ 
suming that all hydrogen is in molecular form and adopt 
the shielding factor given by Draine & Bertoldi (1996), 
their formula 37 (A. 13, see Sec. A.3 for details), and then 
obtain the H 2 abundance from chemical network. The self¬ 
shielding of CO has been calculated in a similar way using 
the shielding functions from Lee et al. (1996). 


3. Results 

3.1. UV radiation and photodissociation rates 

The temperature and density structure of the disk model 
by D’Alessio et al. (1999) is adopted as our fiducial model, 
which was calculated for an accretion rate M = 1O _8 M 0 
yr -1 and viscosity parameter a = 0.01. The disk has a 
radius of 400 AU and a mass of 0.05 M 0 , representative 
of the disk around the T Tauri star LkCal5 for which 
molecular data exist (Qi 2001, van Zadelhoff et al. 2001, 
Thi et al. 2002). Similar sizes or even larger have been 
observed for a number of gas disks, see, e.g., Simon et al. 
(2001): LkCa 15 (650 AU; 435 AU by Qi 2001), DM Tau 
(800 AU), DL Tau (520 AU) and GM Aur (525 AU). These 
are the same disks for which observations of molecules 
other than CO are available and for which chemical models 
have been presented in Paper I. The adopted disk model 
is gravitationally stable according to the Toomre criterion 
up to at least 340 AU (D’Alessio et al. 1999, see also Bell 
et al. 1997 for a discussion on this point). We therefore 
truncated our disk at the radial grid cell in which this 
radius is reached, J? out =400 AU. 

The disk structure is calculated self-consistently for a 
stellar temperature of 4000 K and R *=2 i? 0 (D’Alessio 
et al. 1999). The disk is assumed to be subjected to UV 
radiation from the interstellar radiation field and from the 
star. The Draine (1978) field is used to represent the inter¬ 
stellar field between 912 and 2000 A, extended to longer 
wavelengths as specified by van Dishoeck & Black (1982). 
For the star, initially the same spectral shape is taken, but 
scaled by a factor such that the intensity at R = 100 AU 
incident on the surface of the disk is a factor of 10 4 higher 
than the interstellar field. This choice was motivated by 
Herbig & Goodrich (1986) and was also adopted in Paper 
I. In Sect. 3.4, different assumptions for the stellar radia¬ 
tion field are explored. It is assumed that the additional 
UV radiation does not change the disk structure. The inte¬ 
grated flux in the range 912-4000 A, - the range important 
for the chemistry-, can contribute to a maximum of 10% 
of the total energy for the stars modeled in this paper. 
The UV radiative transfer was calculated using the grain 
properties and extinction curve by Weingartner & Draine 
(2001). 

Figure shows the resulting distribution of the UV 
intensity in the disk for two different wavelengths, with 
overplotted the density and temperature structure. It is 
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Fig. 1. Normalized mean intensities at 1000 and 3400 A for the adopted model. The mean intensities have been nor¬ 
malized to the value in the innermost cell to better visualize the effects of grain scattering and absorption. Overplotted 
are the H 2 number density contours (from surface to midplane: 5xl0 4 , 2xl0 5 , 10 6 , 10', 10 8 , 10 9 ) [cm -3 ] on the left 
and the dust temperature contours (upper left to lower right: 90, 70, 50, 40, 30, 20) [K] on the right. Only one quadrant 
of the 2D flaring disk is shown. 



Fig. 2. Spectrum of the mean intensity J\ at different 
heights in the disk for a radius of R = 105 AU. Numerals 
in the figure indicate the height from the midplane. The 
stellar UV is assumed to be 10 4 times higher than inter¬ 
stellar field at R = 100 AU and have the same shape. 

seen that for longer wavelengths, the radiation penetrates 
deeper into the disk because of the reduced absorption. 
Figure ^ shows the UV spectrum at various heights at 
R = 105 AU. Even at short wavelengths (912 -2000 A), 
the intensity is still comparable to the unattenuated in¬ 
terstellar radiation field down to 35 AU at R ~ 100 AU. 
The longer wavelengths visible photons start to dominate 
the radiation field below 50 AU. The dip at A ~ 2200 A 
is caused by the bump in the extinction curve. 

The effect of radiation on the molecular abundances is 
reflected in the different dissociation or ionization rates. 
For example, Fig. shows the dissociation rates of C 2 H 


(C 2 H -> C 2 + H) and H 2 CO (H 2 CO —> CO + H + H). 
C 2 H is dissociated by lines in between 1080-1530 A and 
H 2 CO is dissociated due to both line (1100 - 1749 A) and 
continuum radiation (912-1505 A). The solid lines show 
the rates calculated from the stellar and interstellar radia¬ 
tion field, obtained via the full 2D radiative transfer. The 
dotted and dashed lines show the rates calculated using 
the depth-dependent dissociation rates of van Dishoeck 
(1988) and NSM, respectively. Van Dishoeck (1988) cal¬ 
culated a plane-parallel PDR model, and obtained the dis¬ 
sociation rates k as functions of depth from the cloud sur¬ 
face. These were fitted by the single exponential decay 
function: 

k = koe~$ A ''. (3) 

The NSM assumes the same function, but with slightly dif¬ 
ferent coefficients. In Aikawa & Herbst (1999) and Paper I, 
the visual extinction A v from the central star is obtained 
at each position in the disk by integrating the density dis¬ 
tribution along the line of sight to the central star, and 
the dissociation rate is obtained with the NSM depth- 
dependent rates. The dissociation rate via the interstellar 
UV field is obtained in a similar way using the depth- 
dependent NSM rates and calculating A v from the disk 
surface. Figure || shows the total dissociation rate via UV 
radiation from the interstellar field and the central star. In 
the dotted and dashed lines, the two components can be 
easily distinguished: radiation from the central star domi¬ 
nates at Z 50 AU, while interstellar radiation dominates 
at Z 50 AU. 

In Fig. [| it is seen that the dissociation rates are signif¬ 
icantly underestimated at Z ~ 50 AU, if the approximate 
depth-dependent rates described above are used. This is 
caused by geometrical effects: the depth-dependent rates 
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Fig. 3. Photodissociation rates of C 2 H and H 2 CO as 
functions of height from the midplane at R = 105 AU. 
Rates are calculated using the radiation field obtained 
via the full 2D Monte Carlo radiative transfer (solid 
lines). For comparison, rates obtained via the approxi¬ 
mate 2D method adopting the depth-dependent function 
of van Dishoeck (1988) (dotted lines ) and the new stan¬ 
dard model (dashed lines) are shown. 

are originally obtained for a plane-parallel cloud, whereas 
the central star irradiates the disk obliquely. For 50 % of 
the photons which are scattered in the disk surface, the 
grazing angle between the disk and the ray path becomes 
larger than the initial value, which causes deeper penetra¬ 
tion into the disk. 

3.2. Vertical distribution of molecules 

Using the photodissociation rates calculated in Sec. 3.1, 
the abundances of the molecules are calculated as a func¬ 
tion of time. Figure [| shows the vertical distribution of 
the molecules at a disk radius of R = 105 AU at 1 10 6 yr. 
The dissociation rates in panels (a) and (b) are obtained 
via the depth-dependent function (Eq. |||) , where the same 
model is used as that presented in Paper I, but plotted at 
a different radius. This model is called the A-2D (approxi¬ 
mate 2D) model hereafter. Two improvements are made in 
this paper: photodissociation rates are calculated via the 
full 2D radiative transfer, and H 2 and CO dissociation via 
stellar UV is included. To disentangle the two effects, two 
models are run. In model ND-2D (no-dissociation, panels c 
and d), the dissociation rates are obtained via the full-2D 
radiative transfer results, but CO and H 2 are not disso¬ 
ciated by the stellar light. Comparison with A-2D (pan¬ 
els a and b) shows that the CO abundances change only 
slightly, caused by the scattering and further penetration 
of interstellar UV photons in full 2D. Most other species 
have lower abundances because the dissociation rate is 
enhanced. The results obtained using the full 2D radia¬ 
tive transfer and H 2 and CO dissociation via the stellar 
UV are presented in panels (e) and (f). CO is dissociated 
even at Z ~ 50 AU in the 2D model, while it is largely 
shielded from UV at Z 60 AU in the other two cases. 
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Fig. 4. Vertical distribution of molecules at R = 105 AU 
and t = 1 10 6 yr. In panels (a) and (b), the approximate- 
2D (A-2D) method is adopted with the depth-dependent 
dissociation rates. The dissociation rates used in panels (c) 
and (d) are determined using the full 2D radiative transfer 
(ND-2D). In both cases the stellar UV is assumed not to 
be energetic enough to dissociate H 2 and CO. In panels 
(e) and (f), dissociation rates are determined using the 
full 2D radiative transfer, and dissociation of H 2 and CO 
via the stellar UV is considered. 


The peak abundances of radicals such as CN and C 2 H (at 
Z ~ 60 AU), are enhanced, because the stellar UV disso¬ 
ciates CO to provide more atomic carbon in the gas-phase 
reaction network. This zone corresponds to the so-called 
“radical zone” found in the chemical networks of PDRs 
(e.g., Jansen et al. 1995a & 1995b, Sternberg & Dalgarno 
1995). The other species benefit from the CO dissociation 
as well; their abundances are higher than in the ND-2D 
model. 

3.3. Column densities 

The column densities obtained at each radius by inte¬ 
grating the vertical distribution of molecules (Fig. |]) are 
presented in Figure [^. The column densities of radicals 
(CN and C 2 H) depend sensitively on the photodissocia¬ 
tion treatment, especially at the inner radii. The column 
densities of CN and C 2 H in the A-2D and ND-2D mod¬ 
els increase with radius because of the lower density and 
lower flux of the destructive stellar UV in the outer re¬ 
gions. Inclusion of CO dissociation provides more carbon 
in the gas-phase reaction network, which significantly en¬ 
hances the abundances of CN and C 2 H. At the inner radii, 
the higher dissociation rates of these radicals are cornpen- 
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sated by a larger carbon supply. It suggests that an accu¬ 
rate calculation of the radiative transfer of the stellar UV 
is important when considering the abundances of radical 
molecules in disks. 



R [AU] R [AU] 


Fig. 5. Radial distribution of column densities at t = 
1 10 6 yr. The dashed lines with open circles refer to the 
A-2D model, the long-dashed lines with the solid triangles 
to the ND-2D model, and the solid lines with closed circles 
are for the full 2D model. 

Photodissociation of CO via the stellar UV at the disk 
surface, however, does not affect the column density of 
CO significantly. It should be noted that Fig. ^ shows 
the molecular abundances relative to hydrogen nuclei, and 
that the absolute number density of CO is highest in the 
layer with Z ~ 30, which contributes predominantly to 
the CO column density. Overall, column densities of other 
species do not significantly depend on the treatment of 
stellar UV. Close to the star, where the effect of the stel¬ 
lar UV field is largest, differences can be seen. The HCN 
abundance close to the star shows that 2D UV transfer is 
important. 

As argued in Paper I and as can be seen in Fig. |], 
molecules are abundant only in layers with certain phys¬ 
ical conditions. Since the mass contained in those lay¬ 
ers does not vary much with radius, column densities 
of molecules such as HCN and H 2 CO show little depen¬ 
dence on radius. The column densities of CO and HCO + 


abruptly change at R ~ 100 AU, inside of which the tem¬ 
perature in the midplane is higher than the sublimation 
temperature of CO (~ 20 K). 

Molecular column densities were also calculated for 
disks with higher and lower accretion rates of 10~ 7 and 
10 -9 Mq yr -1 and correspondingly larger and smaller disk 
masses. For the same reasons stated above, the column 
densities are similar to those in Fig. [|, except for CO and 
HCO + in the inner radius (R 100 AU) and H 2 , whose 
column densities are higher (lower) in more (less) massive 
disks (Paper I). 

3.4. Dependence on the stellar spectrum 



Fig. 6. Comparison of the mean intensity at a radius of 
R = 105 AU for the three stellar spectra: Spectrum A 
(dashed line), Spectrum B (solid line), and Spectrum C 
(dotted line) (see text). 

So far, the UV spectrum from the central star has been 
assumed to be the same as that of the interstellar radia¬ 
tion field, but with an intensity 10 4 times higher than the 
interstellar field at R = 100 AU (hereafter referred to as 
Spectrum A). This spectrum was motivated by Herbig & 
Goodrich (1986) comparing IUE data in the wavelength 
range between 1450 and 1850 A towards different T-Tauri 
stars. In this subsection, the sensitivity of our results to 
the UV spectrum of the central star is investigated by 
calculating two other models with different stellar spec¬ 
tra: the observed spectrum of the T Tauri star TW Hya 
(Spectrum B), and a black body spectrum at T* = 4000 K 
without any excess UV (Spectrum C). Costa et al. (2000) 
fitted the UV continuum from TW Hya, observed with 
the IUE satellite, with a sum of free-free plus free-bound 
emission at 3 x 10 4 K, plus a 4000 K black body spectrum 
for the star and a black body emission at 7900 K cov¬ 
ering 5 % of the stellar surface. We adopt his numerical 
data at the stellar surface as Spectrum B. Note that no 
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observational constraints are available below 1200 A, espe¬ 
cially in the critical 912-1100 A range where H 2 and CO 
are photodissociated. It will be important to obtain such 
data in the future with the FUSE satellite. The third case 
with Spectrum C gives a lower limit to the UV radiation 
from the central star, and may correspond to the situa¬ 
tion for disks around weak-line T Tauri stars (WTTS), 
although the physical structure of such disks would be 
different from those around classical T Tauri stars. As in 
the full 2D model of the previous section, the photodisso¬ 
ciation rates are obtained via full 2D radiation transfer, 
and dissociation of H 2 and CO via stellar UV is taken into 
account. Photodissociation by the interstellar radiation is 
included in all models. Figure [i] shows the spectrum at the 
disk surface at R « 100 AU for the three stellar models. 



Fig. 7. Vertical distribution of molecules at R = 105 AU 
and t = 1 10 6 yr. Panels (a) and (b) are for Spectrum A, 
panels (c) and (d) are for Spectrum B and panels (e) and 
(f) for Spectrum C 


Figures [?] and || show the vertical distribution of 
molecules at R = 105 AU, and the radial distribution of 
column densities for the three models with different stel¬ 
lar UV spectra, respectively. The model using Spectrum B 
shows almost the same distributions and column densities 
as those obtained with Spectrum A. This is not surprising 
since the spectrum of the interstellar radiation field has a 
color temperature of 3x 10 4 K, close to that of the free-free 
and free-bound component in Spectrum B. In the model 
without excess UV (Spectrum C), on the other hand, the 
molecular distributions are far more extended toward the 
surface layers. The lower photodissociation rates at the 
disk surface tend to enhance the column densities, but 
the radical (CN and C 2 H) column densities in the inner re- 



Fig. 8. Radial distribution of molecular column densities 
at t = 1 10 6 yr. Solid lines with closed circles are for 
Spectrum A, dashed lines with open circles for Spectrum 
B, and long-dashed lines with triangles for Spectrum C. 


gions are lowered because of the smaller amount of atomic 
carbon supplied via CO dissociation in the intermediate 
layers. 


4. Molecular line emission 

4.1. Comparison of models 

The comparison of our new results obtained via the full-2D 
radiation transfer with the A-2D model and the ND-2D 
model shows some differences in molecular distributions, 
especially for radicals. 

The question is if these differences are significant with 
the present day observations or whether any differences 
are minimized after solving the radiative transfer in the 
millimeter lines and convolving the data with the telescope 
beam. The CO and HCO + (sub-)millimeter lines are ex¬ 
tremely optically thick in these models. To get more infor¬ 
mation on the emission from deeper layers in the disk, the 
radiative transfer has been run for 13 CO and H 13 CO + as 
well. These species are assumed to follow the abundances 
of CO and HCO + respectively, but their abundances are 
lowered by a factor of 60, similar to the [ 12 C]/[ 13 C] iso¬ 
topic ratio in the solar neighborhood. This overestimates 
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the 13 C0 abundance for the full-2D case since CO is more 
self-shielding compared to 13 CO. 

As in Paper I, the level populations and emission lines 
of different species are calculated using a 2D line radiative 
transfer Monte Carlo code (Hogerheijde & van der Tak 
2000). The calculations assume that the 400 AU radius 
disk is seen at an inclination of 60° at a distance of 150 
pc, appropriate for the case of LkCal5. The results are 
presented in two ways; first, the emission profiles for a 
few characteristic molecules are plotted, and second, the 
ratios of integrated intensities of the lines are compared 
in several tables. In Fig. ^j, lines are shown for CO 3-2, 
HCO+ 4-3, CN 3-2 and HCN 4-3. The CN 3-2 line is 
modeled assuming that the three fine-structure lines are 
well shifted from each other, with each fine-structure line 
a combination of three unresolved hyperfine components. 
Therefore, the CN 3-2 line shows a regular double-peaked 
profile instead of a combination of three double-peaked 
structures summed at slightly shifted velocities. 


Table 1. The ratios of the integrated intensities for the 
full-2D, ND-2D and A-2D models (see text). The line 
emissions have been convolved with a beam equal to the 
size of the disk on the sky (5.33"). 


Molecule 

line 

v [GHz] 

full-2D 

A-2D 

ND-2D 

A-2D 

CO 

6-5 

691.473 

0.88 

0.97 

CO 

3-2 

345.796 

0.92 

0.98 

CO 

2-1 

230.538 

0.93 

0.98 

13 CO 

3-2 

330.587 

0.98 

1.03 

13 CO 

1-0 

220.399 

1.05 

1.08 

HCO+ 

4-3 

356.734 

0.94 

0.93 

HCO+ 

1-0 

89.189 

0.91 

0.89 

h 13 co+ 

4-3 

346.999 

0.94 

0.93 

h 13 co+ 

1-0 

86.754 

0.98 

0.96 

CN 

3-2 

340.248 

10.1 

0.75 

CN 

2-1 

226.875 

5.81 

0.85 

HCN 

4-3 

354.506 

1.21 

0.81 

HCN 

1-0 

88.632 

1.19 

0.85 


6.0 



0.0 


1.5 





0.0 



Fig. 9. Comparison of emission lines for the full-2D radia¬ 
tive transfer model (solid lines), the 2D radiative trans¬ 
fer model without CO and H 2 dissociation (ND-2D; dot- 
dashed) and the approximate 2D model (A-2D; dashed 
lines). The CO 3-2, HCO+ 4-3, CN 3-2 and HCN 4-3 
lines are convolved with a beam equal to the size of the 
source on the sky (5.33 "). 


In Table |l], the ratios of the integrated intensities of the 
full-2D and A-2D as well as the ND-2D over A-2D models 
are given for a few lines. All models refer to a beam size 
equal to the apparent size of the disk at a distance of 
150 pc. In both Fig. [ 9 ] and Table |l| the CO emission is 
slightly lowered in the full-2D case due to the dissociation 
by stellar UV. 


Due to the dissociation of CO leading to atomic car¬ 
bon, the CN abundance and emission is enhanced up to 
an order of magnitude. This is clearly seen in Fig. ^j, where 
the A-2D and ND-2D models have similar emission pro¬ 
files, but full-2D is much higher. HCN is linked to CN and 
its abundance is enhanced as well; however, in contrast 
with CN, which can only be dissociated at wavelengths 
less than 1000 A, HCN can be easily destroyed again by 
photodissociation, resulting in both the abundances and 
line emission to rise only mildly. Among the three models, 
the ND-2D model gives the weakest HCN emission. 

4.2. Dependence on stellar radiation field 

In the previous section it was shown that the chemical 
abundances need a good description of the mean intensity 
of the UV radiation field. The column densities of radicals 
such as CN are especially sensitive to the UV spectrum. 
The models with excess UV radiation (Spectrum A and 
B) have a flat column density distribution throughout the 
disk, while the model without excess UV (Spectrum C) 
has an increasing CN column density toward the outer ra¬ 
dius. Figure [TTi| shows line profiles of HCN, HCO + and CN 
for the three different stellar spectra. The lines of HCO + 
nicely reflect the column densities shown in Fig. ||; the 
model with Spectrum C gives the strongest emission. The 
HCN column densities have similar values at larger radii, 
which is reflected in the 4-3 emission lines. At smaller radii 
the column of the model with Spectrum C dominates by a 
factor of 5. This is reflected in the higher emission in the 
wings of the line profile. The excess emission at |v| ~4 km 
s” 1 is due to the undersampled grid in the inner region, 
but this does not effect the integrated emission. The CN 
emission is discussed in more detail in Sec. |T~^ . 

For a broader overview of the calculated emission lines, 
the integrated intensities are given in Table || for several 
species and lines. These results were again calculated for 














10 


van Zadelhoff et al.: Chemical effects due to 2D UV transfer 





v (km s ^ v(kms 1 ) v (km s ) 


Fig. 10. Comparison of HCO + 4-3 (left), HCN 4-3 (middle) and CN 3-2 (right) emission lines for the three stellar- 
spectra. The dash-dotted lines are for Spectrum A, the solid lines for Spectrum B, and the dashed lines for Spectrum 
C. The spectra have been convolved with a beam of 5.33" 


Table 2. Integrated intensities for the three stellar spec¬ 
tra (A, B and C) convolved with a beam of 5.33". 


molecular line 

A 

[K km s" 1 ] 

B 

[K km s -1 ] 

C 

[I< km s -1 ] 

CO 

6-5 

11.19 

11.79 

12.65 

CO 

3-2 

16.92 

17.62 

18.74 

CO 

2-1 

17.86 

18.52 

19.58 

13 CO 

3-2 

7.93 

8.05 

8.28 

13 CO 

2-1 

7.84 

7.93 

8.09 

13 CO 

1-0 

5.28 

5.31 

5.38 

HCO+ 

4-3 

3.38 

3.72 

4.45 

HCO+ 

1-0 

2.07 

2.31 

2.93 

h 13 co+ 

4-3 

0.10 

0.11 

0.14 

h 13 co+ 

1-0 

0.04 

0.05 

0.06 

CN 

3-2 

0.68 

0.66 

0.15 

CN 

2-1 

1.94 

2.00 

1.01 

HCN 

4-3 

0.35 

0.37 

0.46 

HCN 

1-0 

0.56 

0.63 

0.77 


the source seen with a beam of 5.33". The CO emission 
lines show only a small dependence on the stellar UV field 
(~ 10%), where the case without excess UV (Specrum C) 
shows the highest integrated flux. The 13 CO lines show an 
even smaller difference since most of the emission comes 
from the warm intermediate layers and is not affected 
by dissociation in the surface layer. HCO + and its main 
isotoponrer H 13 CO + have a similar behavior, but with 
a larger difference (~ 30-50%). The HCO + emission is 
strongest for the star with no excess UV (Spectrum C), 
and is slightly stronger for Spectrum B than A. 

4.3. NLTE effects on the CN line emission. 

Special attention should be paid to the CN lines, because 
they are among the strongest lines observed in several 
disks, and because LTE excitation is not a good approxi¬ 
mation for this molecule as shown below. The integrated 
intensities of the CN lines are larger by a factor of 2 to 3 


(Table ||) in the models with excess UV (Spectrum A and 
B), even though the CN column densities are a factor of 
two higher in the outer regions of the disk for the model 
without excess UV (Spectrum C) (Fig. ||). This appears 
counterintuitive since the outer disk, due to its larger sur¬ 
face, should dominate the line emission. This is not the 
case since most of the additional CN in the model with 
Spectrum C is located in the upper most layers of the 
disk where the density is so low that the levels become 
sub-thermally excited. To show the effect of the NLTE 
excitation, the profiles of the CN 3 — 2 line for the three 
stellar spectra are plotted for NLTE and LTE populations 
in Fig. [H]. In the LTE case the emission from the model 
with Spectrum C dominates the lines from Spectrum A 
and B as expected when considering the column densities 
(Fig. ||). For the NLTE calculations, the order of the emis¬ 
sion reverses since most of the CN column in the model 
with Spectrum C is at low density. 



v (km s ) v (km s ) 

Fig. 11. Comparison of the CN 3 — 2 emission line 
for NLTE (left) and LTE (right) populations. The dash- 
dotted lines are for Spectrum A, the solid lines for 
Spectrum B, and the dashed lines for Spectrum C. The 
spectra have been convolved with a beam of 5.33" 

In addition, a few of the hyperfine lines of the J = 
1 — 0 transition experience small masing effects due to a 
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population inversion. This makes NLTE excitation effects 
important for all CN levels in these disk environments. 
The modeling and interpretation of any of the CN levels 
should therefore be done with a detailed radiative transfer 
study. 

4.4. High resolution simulation of CN/HCN 

The CN and HCN column densities (Fig. ||) show clear 
differences for the different radiation fields. Their abun¬ 
dance ratio in the disk is a sensitive function of radius 
due to the different radial gradients of the column den¬ 
sities of CN and HCN. In the future, telescopes like the 
Atacama Large Millimeter Array (ALMA) will be able 
to probe this ratio spatially. In order to simulate ALMA 
observations, the molecular line intensities are calculated 
assuming a FWHM gaussian beam of 0.3", which is cho¬ 
sen based on the spatial resolution of the model grid and 
which will easily be reached with ALMA. In Fig. |l2| c the 
velocity integrated intensity ratio along the major axis of 
the disk is plotted for spectra B and C. The intensity ratio 
for spectrum B peaks close to the star and drops to a 15 
times lower ratio at large radii. Spectrum C shows the op¬ 
posite behavior rising by a factor 4 from close to the star to 
larger radii. The emissions at these large radii are similar 
(Fig. |l2|b,c), so that high resolution is essential for distin¬ 
guishing the effects of different UV radiation fields. For 
a better understanding of the emission profile, the beam 
averaged column densities and integrated intensities for 
CN and HCN are given in Fig |l^ a and 6, respectively. 
The drop of intensities at ~ 2.5" is due to the partial 
filling of the beam. In both cases the ratio of the inten¬ 
sities for spectrum B and C is smaller than the ratio of 
the column densities by roughly an order of magnitude; 
in the inner most region with spectrum C, for example, 
the column density ratio of CN to HCN is almost unity , 
while the intensity ratio is <; 0.1. The expected sensitivity 
of ALMA at 345 GHz in a 0.3" beam is ~0.3 K km s -1 
in 1 hr. Thus, the HCN and CN distribution in the inner 
0.6” can readily be imaged and the Spectrum B and C 
models should be easily distinguishable. Imaging of CN in 
the outer part of the disk will require longer integration 
times. 


5. Discussion 

The calculated abundances and column densities show or¬ 
ders of magnitude variations between different models. 
The inclusion of ND-2D UV radiative transfer instead 
of an approximate 2D transfer decreases the abundances 
of species like HCN, CN and CS. When CO dissociation 
via stellar UV is taken into account, a large amount of 
atomic carbon is produced, which is converted to radicals 
(CN, C 2 H) and direct reaction products (CS, HCN). The 
adopted stellar spectrum has a clear effect on the column 
densities of most species except for CO. Especially the 
CN and HCN column densities show different slopes with 



0 12 3 


distance along major axis (arcsec) 


Fig. 12. a) The beam averaged column densities of CN 
and HCN. b) Integrated intensities for the CN 3-2 and 
HCN 4-3 transition, c) Ratio of the integrated intensity 
of CN 3-2/HCN 4-3. The modeled intensities and columns 
were calculated using a beam of 0.3". The solid lines rep¬ 
resents the emission of a star with stellar UV (Spectrum 
B), the dashed lines a T-Tauri star with no excess UV 
(Spectrum C). The thick lines in the plots (a) and (b) 
represent the CN and the thin lines the HCN column- 
densities and intensities. 


increasing radius. These variations should be detectable 


with sufficient spatial resolution as modeled in Fig. |12 


The effect of UV radiative transfer and the stellar UV 
spectrum are less clear in the integrated line intensities 
with larger beams for the following reasons. First, the tele¬ 
scope beams available, now and in the near future, average 
the disk over the entire beam. Second, most lines observed 
to date are the brightest lines available. A large number of 
these are optically thick; CO for instance is only capable 
of probing the upper layers of the disk. Third, some lines 
suffer from sub-thermal excitation. The NLTE effects can 
be severe in the optically thin upper layers of the disk, 
while the LTE intermediate layers are at a large optical 
depth. 

A comparison between different millimeter line radia¬ 
tive transfer codes (van Zadelhoff et al 2002) shows that 
the standard deviation on the level populations can reach 
12% for the higher levels in NLTE. These calculations were 
performed for an extreme combination of optical depth 
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(t) and NLTE. Even though the r in disks is very high, 
the levels are not as far from LTE as in the above men¬ 
tioned calculation. It is therefore safe to assume that the 
molecular line calculation in this paper (Fig. 9 and 10) has 
‘error-bars’ similar or less than 12%. This is smaller than 
the error-bars on the observations, which are estimated to 
be roughly 20%. 

In Paper I, the model results underestimated the ob¬ 
served molecular abundances of species such as CN and 
HCN, and overestimates CO in most disks. One of the 
motivations of the present study is to check if the inclu¬ 
sion of CO dissociation via stellar UV could resolve these 
discrepancies. In comparison to the results presented in 
Paper I, both the CN emission and the CN/HCN ratio 
has significantly increased using the full-2D UV radiative 
transfer code. The column density of C 2 H is also increased. 
However, the CO column density itself is not much af¬ 
fected by dissociation via the stellar UV. The column den¬ 
sities of non-radical species such as HCN are almost the 
same as in Paper I and are thus lower than observed in 
some sources. 

There are several important physical processes which 
need to be investigated in future studies (see also §2.1). 
The UV transfer was calculated here using a particular 
set of grain scattering properties. The mean intensities 
derived in the disk change when different properties are 
adopted. Since the abundances do not depend linearly on 
the mean intensity, they are not expected to differ much 
as long as the gas/dust ratio and grains sizes are simi¬ 
lar to those in interstellar clouds. If the dust sizes have 
increased to mm-sized grains, this would alter both the 
dust properties and the dust settling. This may allow more 
UV photons to reach deeper into the disk, enhancing the 
dissociation of CO. 

The inclusion of X-rays could significantly change the 
chemical abundances (Aikawa & Herbst 1999; 2001). The 
secondary UV field generated by the excited H 2 could re¬ 
sult in a higher abundance of radicals at lower disk heights. 
This would enhance the total line emission due to both 
the higher CN and HCN column densities and the in¬ 
creased excitation of the molecules at higher densities. 
Hydrodynamics such as accretion and turbulent mixing 
would smooth out molecular stratification to some extent, 
and transport of radical species would affect the chemical 
processes in neighboring regions. On the observational side 
spatially and spectrally resolved observations are needed 
to find unique tracers of different processes. 

6. Conclusion 

The main conclusions from this work are: 

— The transfer of stellar ultraviolet radiation has to be 
treated in a two-dimensional manner for a correct 
abundance determination of radicals like CN and C 2 H. 

— The photodissociation of CO by the stellar radiation 
supplies additional atomic carbon, which leads to sig¬ 
nificantly higher abundances of molecules, especially 


of radicals like CN and C 2 H. CO itself shows no de¬ 
tectable dependence on the stellar radiation field, mak¬ 
ing it a poor tracer of the chemical processes in disks. 

— The abundances of the molecules in the ‘radical re¬ 
gion’ are also sensitive to the adopted stellar UV spec¬ 
trum. This effect will not necessarily be reflected in 
the current observed line intensities, however, since the 
lines may be sub-thermally excited in the surface lay¬ 
ers where the abundances are most enhanced. Spatially 
resolved observations by ALMA are needed to distin¬ 
guish the different scenarios. 

— Observational constraints on the fluxes from T Tauri 
and Herbig Ae stars at ultraviolet wavelengths, in par¬ 
ticular in the critical 912-1100 A region, are urgently 
needed. 
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Appendix A: Axi-symmetric 2D continuum 
radiative transfer 

A two-dimensional (2D) axi-symmetric Monte Carlo code 
was set up to calculate the transfer of ultraviolet radiation 
by absorption and scattering due to grains. In the calcu¬ 
lation, the photons are actually followed in 3D, but by 
applying axial symmetry the grid can be rotated at any 
time in 4> such that the photon comes back to the <F=0 
plane. This causes the photons to follow hyperbolic paths 
in the 2D <f>=0 plane . The main reason for adopting a 
2D rather than 3D approach is the increase in accuracy in 
the solution for the same number of photons or amount of 
computer time. 


A.l. Photon propagation 


Each photon is started at the source of radiation with an 
intensity given by the boundary conditions for the source. 
These boundary conditions are discussed in Sec. A.4. A 
photon package can undergo two different processes: ab¬ 
sorption and scattering by dust. A model photon repre¬ 
sents the amount of energy in a wavelength bin, given 
by the total surface of the source multiplied by the flux 
through this surface divided by the number of photons in 
this wavelength bin. Per wavelength bin, one needs suffi¬ 
cient resolution to sample all possible angles of emission. 
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In general the total number of photons per wavelength bin 
6X is given by 


N P hot,s\ = Cl x n(j> x n6 x nR x nz, 


(A.l) 


where nR and nz are the number of cells in the R— and 
z—directions, respectively, and ncj) and nd are the direc¬ 
tional resolutions. One can multiply this by an arbitrary 
constant (Cl > 1) to improve on the statistics. For a rea¬ 
sonably sized grid (nR = nz = 70), 10 6 photons suffice 
in practice per wavelength bin. In general 34 wavelength 
bins are used in between 912 A and l pm. 

The initial energy stored in each photon package is 
given by 


A(0, A) 


F x xS 

Nphot,S\ 


(A.2) 


where F\ is the flux entering the system and S the total 
surface it passes through. 

At the surface of emission source, the scattering optical 
depth is first calculated using a random number generator 


Tscat = -ln(l - C) 


(A.3) 


with C, a random number between 0 and 1. This optical 
depth can be converted to a total absorption optical depth 


Tabs — F \ X Tscat-, (A. 4) 

where I\x is the fraction of absorption compared to scat¬ 
tering. The photon package travels through the grid until 
it reaches the optical depth ( r a bs ) at which it scatters. In 
the mean time, its intensity drops according to 


Ii(s + As, A) = Ii(s, A)e Arx , 


(A.5) 


where Ary is related to the physical properties through 


Ata 


1 

2.51ogi 0 (e) 
gd(R,z) 
2.51ogio(e) 


x AA a 

As x n(H + 2H 2 ) 

X --- 7 T - 

1.59-10 21 


Aa 

A v 


(A.6) 

(A.7) 


with Aa the visual extinction at wavelength A, gd(R,z) the 
gas-to-dust ratio compared to 100, As the path length 
traveled within a grid cell and n(H + 2 H 2 ) the total num¬ 
ber density in the grid cell. The relation n(H)+ 2n(H 2 )= 
1.59-10 21 A„ is taken from Savage et al. (1977). When the 
photon package reaches the point of scattering, the en¬ 
tire package changes direction according to the Henyey- 
Greenstein phase function (Henyey-Greenstein 1941, Witt 
1977): 


P(coso,5a) 


_ 1 ~9x _ 

47 r[l + gf — 2gx cos a] 3 / 2 


(A.8) 


with gx =< cosa >, the mean scattering angle. The an¬ 
gles under which the photon scatters are in reference to 
the original direction of motion [defined with an ( x , y, z) 


coordinate system, such that z is in the direction of the 
photon propagation prior to scattering] 


a = arccos 


(3 = tt(2C - 1) 


(1 + 5a) _ 


1 -9l 


1 - 5a + 2<?aC 



with a the angle between the z-axis and the new direction 
(z n e w) and (3 the rotational angle in the x — y plane (Fig. 
A.l), with £ a random number between 0 and 1. 


Z 



Fig. A.l. Scattering of a photon propagating in the direction 
(4>,6) through angles a and /3. The new photon direction will 
lie on a cone defined by a and f3. 


These angles are transformed back to the rest frame 
(. X,Y,Z). After the new direction is calculated the next 
step is determined using equations A.3 and A.4. The pho¬ 
ton continues its path until it reaches the outer boundary 
of the computational grid or has a negligible intensity de¬ 
fined by the user. 


A.2. Computation of the mean intensity 

The mean intensity is determined in each cell through 
summation of all photons passing through the grid cell. 
As each photon represents an amount of energy which 
diminishes as it traverses a (dense) medium, the energy 
added to the cell is the mean energy per second of the 
photon along the path As traveling with the velocity of 
light. The mean energy density ux(iR,iz) of the radia¬ 
tion field is then calculated by dividing the total summed 
energy by the volume V of the gridcell iR, iz. 


ux(iR , iz) 


1 \ ' j As (1 - e At ) 

V(iR,iz)^—' 1 c At 


(ATI) 
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The mean intensity J\(iR,iz) in the grid cell is directly 
related to the mean energy density u\(iR,iz) according 
to 

J\(iR,iz) = —u\[iR,iz). (A.12) 

47T 

The mean intensity is the only output needed in this 
application, so the above description of scattering is suf¬ 
ficient. If the code is to be applied for both intensity and 
polarization, the four Stokes parameters have to be cal¬ 
culated (e.g., Whitney 1991). This is in principle readily 
implemented in the code, but has not been done here. 


A.3. H ‘2 and CO self-shielding 


Photo dissociation of H 2 proceeds through absorptions in 
the Lyman and Werner bands in the 912-1100 A region. 
About 10-15 % of all excitations lead to the dissociation 
of H 2 . Since this process proceeds through line excitation 
the absorption lines become optically thick once the H 2 
column in each level reaches ~10 14 cm' 2 . As this is eas¬ 
ily reached for the objects of interest in this paper, self¬ 
shielding of H 2 has to be taken into account. Only a lim¬ 
ited amount of continuum radiation can be absorbed by 
H 2 , leaving a large flux available for other species to be 
dissociated. The equation used to calculate the shielding 
due to H 2 is taken from Draine & Bertoldi (1996), then- 
equation 37 is shown here as A.13 . The equation describes 
the self-shielding of the H 2 molecule through the Lyman 
and Werner series: 


Among the other molecules, only CO is affected sig¬ 
nificantly by self-shielding. In the other cases the dust 
extinction combined with H 2 extinction, if the dissocia¬ 
tion is at similar wavelengths (e.g., CN photodissociation 
or C photoionization), is sufficient to describe the rates of 
ionization and dissociation. CO shielding is treated in a 
similar way as H 2 with the difference that instead of using 
an analytical solution, the values are taken from a table 
by Lee et al. (1996), in which shielding factor is given as 
a function of H 2 column density, CO column density and 
A w . Calculation of the shielding factor thus require the CO 
abundance in each grid cell. In principle, iteration on the 
chemistry and radiative transfer is needed to calculate the 
precise CO dissociation rates. In the case of circumstellar 
disks, a large column is reached for each single cell due to 
the extremely high total columns, especially in the inner 
parts of disks. This causes the stellar dissociation rates to 
depend only on the mean value of the CO abundance in 
each cell. A constant CO abundance of CO/H 2 =7 • 10 -5 
is assumed when calculating the CO self-shielding. If one 
is interested in a medium for which the precise positional 
change in the CO shielding is important, the spatial grid 
near the boundary at which the flux enters has to be in¬ 
creased and the chemistry and radiative transfer have to 
be iterated. Only a small chemical network is needed for 
this purpose, since the CO abundance depends mostly on 
a few simple molecules and atoms. 

A.4. Sources of radiation 


fs( H 2 ) 


0.965 

(1 — x/b5) 2 


0.035 

y/l — X 


x e 


-8.5x10" 


(A.13) 


with x= ^ and b ,5 = b ^ and f s the shield¬ 

ing of the UV in the H 2 lines. Since H 2 dissociates through 
lines, only a fraction of the continuum in between 912 - 
1110 A is absorbed. 

This equation is valid for each photon package j trav¬ 
eling through a gridcell since the formulae is defined as the 
drop in intensity of the radiation held due to the column 
of H 2 . First the mean value of the function f a , for each 
package during its step is calculated, where the mean is 
taken over the column density N of the step: 


F(j,iR,iz) 


Sm fs(j,iR,iz)dN 

Ini 


(A.14) 


This then leads to a mean value for the grid cell 
through: 


F ce n(iR,iz) 


E np(iR,iz) 
3 =0 


F(j,iR,iz)As(j,iR,iz) 


E np(iR,iz) 
3=0 


As(j,iR,iz ) 


(A.15) 


with np(iR , iz) the number of photon steps taken in the 
grid cell defined by iR and iz and As the length of the 
step of the photon. The value of F ce u is then tabulated for 
each cell and applied to the intensities between 912-1100 
A, which suffered dust extinction only. 


Radiation can come from several sources. Discussed here 
are two sources of radiation: interstellar radiation and stel¬ 
lar radiation. Since all boundary conditions consist of a 
surface through which flux passes, there is a general de¬ 
scription for the angle distribution of the photon packages. 
These angles, 

0 = cos"Vl — C (A. 16 ) 

<F = 2ttC (A.17) 


are given in the frame in which the z-axis lies perpendic- 
ularto the t angen t plane of the source surface. Equations 
( A.16| ) and ( A.17 ) need to be transformed to calculate the 
angles in the plane where its z-axis is parallel to the Z- 
axis of the problem (Fig. |A.2| ). The equations have been 
taken from Yang et al. (1995) and describe the isotropic 
emission from a surface. 


— Interstellar field: The interstellar field is assumed to be 
isotropic over 47 t steradian while the computational 
domain is a rectangle in R — z space (Fig. |A.2| ). To 
calculate the mean intensity, a virtual sphere with a 
radius equal to that of the upper right corner is taken 
in order to represent the angle distribution of the in¬ 
coming Held. The starting positions of the photons are 
randomly taken at the surface of the sphere with a 


direction given by Equations A.16 and A.17. These 
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R 


Fig. A.2. The grid represents the computational domain with 
a virtual shell around it. The dashed line denotes the the frame 
to describe the angle distribution of incoming photons from 
the boundary. The solid curved line is the path of the photon- 
package, the dot-dashed line is its virtual path. 


directions are in the frame denoted by z. For the cal¬ 
culation, the directions are transformed to the com¬ 
putational frame denoted by Z. The photon packages 
are placed at the edge of the computational rectangle 
after w hich the photons follow the routine described 
in Sec. A.l. At each grid cell ( iR,iz ) through which 


the photon passes, the mean intensity is updated. In 
general the Draine interstellar radiation field with the 
van Dishoeck & Black (1982) extension for A > 2000 
A is used which can be scaled by a factor if required. 
Stellar light: Since the medium in our problems is gen¬ 
erally much larger than the stellar radius, the star can 
safely be assumed to be a point source. However, the 
total energy emitted does depend on the radius of the 
star. In this paper we assumed R+=2Rq The far-field 
expression for the mean intensity in a medium with 
negligible absorption is, 


M r ) = 


h (R 


(A.18) 


where I\ is the intensity at the surface of the star. It 
is a good approximation in the region of R j> 5i?*, 
which we are concerned with. For problems where the 
mean intensities close to the star are important, emis¬ 
sion from the real stellar surface should be taken into 
account. 
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